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Abstract  In  this  work  we  consider  the  Navier-Stokes-Korteweg  equations,  a  diffuse- 
interface  model  describing  liquid- vapor  phase  transitions.  A  numerical  scheme  for 
this  model  is  constructed  based  on  functional  entropy  variables  and  a  new  time  inte¬ 
gration  concept.  The  fully  discrete  scheme  is  unconditionally  stable  in  entropy  and 
second-order  time-accurate.  Isogeometric  analysis  is  utilized  for  spatial  discretiza¬ 
tion.  The  boiling  problem  is  numerically  investigated  by  making  proper  assump¬ 
tions  on  transport  parameters  and  boundary  conditions.  Compared  with  traditional 
multiphase  solvers,  the  dependence  on  empirical  data  is  significantly  reduced,  and 
this  modeling  approach  provides  a  unified  predictive  tool  for  both  nucleate  and  film 
boiling.  Both  two-  and  three-dimensional  simulation  results  are  provided. 


1  Introduction 


Boiling  is  a  thermally  induced  phase  transition  process  in  which  new  liquid-vapor 
interfaces  are  generated  in  a  bulk  liquid  region  [2] .  It  is  an  extremely  effective  mech¬ 
anism  in  energy  transfer  and  is  widely  used  in  energy  conversion  facilities.  Despite 
its  importance  in  industry,  the  fundamental  mechanism  of  boiling  is  still  not  well 
understood  [2] .  A  predictive  model  for  boiling  is  highly  desired  for  engineering  de¬ 
signs.  Film  boiling  is  regarded  as  most  amenable  to  modeling,  since  its  governing 
mechanism  is  principally  the  Rayleigh-Taylor  instability.  However,  existing  simula¬ 
tions  all  start  with  a  preexisting  perturbed  flat  interface  as  the  initial  condition  [8]. 
In  other  words,  none  of  those  methods  captured  the  film  generation  process.  On  the 
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other  side,  very  few  simulations  of  nucleate  boiling  have  been  performed  because 
more  physical  mechanisms  are  involved  in  this  phenomenon. 

Traditional  interface-tracking  and  interface-capturing  methods  are  designed  based 
on  geometrical  information  of  existing  interfaces.  This  is  perhaps  the  reason  why 
these  methods  become  intractable  for  phase  transition  phenomena.  Phase-field 
or  diffuse-interface  methods  were  proposed  as  an  alternative  interface-capturing 
method,  that  use  thermodynamic  state  variables  to  distinguish  different  phases  [1]. 
The  solid  mathematical  and  thermodynamic  foundations  of  phase-field  models  al¬ 
low  them  to  describe  these  complicated  phenomena  without  resorting  to  model¬ 
ing  “tricks.”  The  initial  instantiation  of  phase-field  methods  is  the  Navier-Stokes- 
Korteweg  equations,  which  are  constructed  based  on  the  van  der  Waals  theory  [1,4]. 
In  the  past  decades,  this  theory  has  been  developed  further  [3],  and  a  rational  ther¬ 
momechanical  framework  for  the  Navier-Stokes-Korteweg  equations  has  been  pre¬ 
sented  very  recently  [10]. 

For  phase-field  problems,  the  non-convexity  of  the  entropy  function  precludes 
the  possibility  of  directly  applying  many  existing  robust  numerical  methodologies 
[11].  To  overcome  the  challenges  posed  by  the  non-convexity  of  the  entropy,  first, 
functional  entropy- variables  are  introduced  to  construct  an  entropy- stable  spatial 
discretization  [9,  10].  Second,  to  develop  a  stable  temporal  scheme,  we  adopt  the 
methodology  based  on  special  quadrature  rules  [5,  9].  This  time  integration  concept 
can  be  viewed  as  a  second-order  modification  to  the  mid-point  rule.  The  modifica¬ 
tions  are  designed  so  that  the  temporal  approximation  is  provably  entropy  dissipa¬ 
tive.  Since  this  temporal  scheme  does  not  require  convexity  for  the  entropy  function, 
it  is  anticipated  to  be  applicable  to  many  more  general  problems. 


2  The  Navier-Stokes-Korteweg  equations 

We  consider  a  fixed,  connected,  and  bounded  domain  Q  C  M3 .  The  time  interval  of 
interest  is  denoted  (0,T),  with  T  >  0.  The  dimensionless  Navier-Stokes-Korteweg 
equations  are  considered  in  the  space-time  domain  Q  x  (0,  T)  as 

^+V-(pu)=0,  (1) 

+  V-  (pu0u)  +  Vp-  V-  T  —  V-  g  =  pb,  (2) 

dt 

)  +  V  •  ((pE  +  p) u  -  (t  +  S)u)  +  V  q  +  V  II  =  pb  u  +  pr.  (3) 
dt 

In  the  above  equations,  p  is  the  density,  u  is  the  velocity,  E  is  the  total  energy,  p  is 
the  thermodynamic  pressure,  T  is  the  viscous  stress,  g  is  the  Korteweg  stress,  q  is 
the  heat  flux;  JI  is  the  interstitial  working  flux  [3,  10],  b  is  the  prescribed  body  force 
per  unit  mass,  and  r  is  the  heat  source  per  unit  mass.  The  constitutive  relations  are 
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wherein  0  is  the  temperature,  K  is  the  conductivity,  Re  is  the  Reynolds  number,  We 
is  the  Weber  number,  y  is  the  heat  capacity  ratio,  and  t  is  the  internal  energy  density 
per  unit  mass.  The  mathematical  entropy  function  H  and  the  local  Helmholtz  free 
energy  Wioc  are  defined  as 
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In  three  dimensions,  the  vector  of  conservation  variables  is 

Ur  =  [Ui,U2,U3,U4,U5]  :=  [p,pui,pu2,pU3,pE]. 


Due  to  the  appearance  of  the  gradient- squared  term,  the  mathematical  entropy  func¬ 
tion  H  is  no  longer  just  an  algebraic  function  of  the  conservation  variables,  but 
rather  it  is  a  functional  of  the  conservation  variables.  We  define  the  entropy  vari¬ 
ables  Vr  =  [Vj ,  V2,  V3 ,  V4,  V5]  as  the  functional  derivatives  of  H  with  respect  to  U: 

SH 

Vi[Svi]  = -j^y[Svi],  i=l,. ..,5, 

wherein  8\T  =  [5vi ,  8v 2, 8v 3, 5v 4, <H>s]  are  the  test  functions.  The  entropy  variables 
V  can  be  written  explicitly  as 

V1[5v1l  =  i(^-!^)«v1  +  T‘Vp.vSv1, 

Vj[8vi\  m  1  8vi,  i  =  2,3,4,  y5[5v5]  = -l5v5, 

wherein 

Vloc  =  -2P  +  ^01n  (,!p)  -  27(7TT)  0  (ln(e)  “  1}  +  27(1 —p) 

is  the  local  electrochemical  potential.  Inspired  from  the  form  of  Vj ,  we  introduce  a 
new  independent  variable  V  as 
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The  fundamental  thermodynamic  relation  between  p  and  Vioc  allows  us  to  express 
p  in  terms  of  V  as 


p  =  pV0-p'Flo 


p\u[ 


We 


p0V- 


(4) 


Making  use  of  the  relation  (4),  the  original  strong-form  problem  (l)-(3)  can  be 
rewritten  as 
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The  new  strong-form  problem  (5)-(8)  is  an  equivalent  statement  of  the  original 
Navier-Stokes-Korteweg  equations  (l)-(3). 


3  The  fully  discrete  scheme 


The  time  interval  (0,  T)  is  divided  into  Nts  subintervals  (tn,tn+\)9  n  =  0,  •  •  •  ,Nts  —  1, 
of  size  Atn  =  tn+\  —  tn.  We  use  the  notation 
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to  represent  the  fully  discrete  solutions  at  the  time  level  n.  We  define  the  jump  of 
density,  linear  momentum,  and  total  energy  over  each  time  step  as 

IPnl  :=  Pn+ 1  -  Pn  .  —  Pn+lU«+l  “  PnU«- 

[pnE(Pn,»-l  0hn)}  ■  =  (pfU(P*+ 1  ,  O  ~  (P*U(P*+1  , 

+(p^)(pB\i,0*+. )  -  (p**U(p£,  e*+i ) 
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With  the  jump  operators  defined  above,  the  fully  discrete  scheme  can  be  stated  as 
follows.  In  each  time  step,  given  and  the  time  step  Atn,  find  Y^+1  such  that  for 
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In  our  work,  Non-Uniform  Rational  B-Splines  (NURBS)  basis  functions  are  used  to 
define  ^  as  well  as  the  computational  domain.  Consequently,  this  approach  may 
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be  considered  as  isogeometric  analysis  method  [7].  The  main  results  of  the  fully 
discrete  scheme  (9)-(12)  are  stated  in  the  following  two  theorems. 

Theorem  1.  The  solutions  of  the  fully  discrete  scheme  (9)-(12)  satisfy 


■L 
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dVx  <  0. 


Theorem  2.  The  local  truncation  error  in  time  &  (t)  =  \Op  (t) ;  0„  (t) ;  @e  (t)  j  can 

be  bounded  by  \0(tn)  |  <  KAt% I5  for  all  tn  G  [0,  T],  where  K  is  a  constant  indepen¬ 
dent  of  Atn  and  I5  =  (1;  1;  1;  1;  \  )T . 


The  proofs  of  the  above  two  theorems  can  be  found  in  [10].  Theorem  1  states  that  the 
method  is  unconditionally  entropy  stable,  because  d3Vioc/ dp3  >  0  and  d^H/dO3,  < 
0,  which  follow  from  properties  of  the  van  der  Waals  fluid.  Theorem  2  establishes 
the  second-order  time-accuracy  of  the  method. 


4  Boiling 


To  obtain  successful  boiling  simulations,  there  are  several  additional  modeling  con¬ 
siderations.  First,  the  transport  parameters  need  to  be  density  dependent  in  order  to 
differentiate  the  properties  of  the  liquid  and  vapor  phases.  In  our  simulations,  the 
dimensionless  viscosity  coefficient  and  the  dimensionless  conductivity  are  modeled 
as 


a  =  cfyp,  K=cbK°iip , 

with  C^°l1  and  C^otl  being  constants  independent  of  p.  Second,  the  gravity  ef¬ 
fect  need  to  be  taken  into  account  to  generate  buoyancy.  The  dimensionless  body 
force  b  is  chosen  as  b  =  (0;0;  — 0.025)r  for  the  three-dimensional  case  and  b  = 
(0;  — 0.025)r  for  the  two-dimensional  case.  Third,  the  ninety-degree  contact  angle 
boundary  condition  is  used  for  the  density  variable,  and  the  slip  boundary  condition 
is  applied  to  the  velocity.  To  specify  the  boundary  condition  for  T5  =  — 1/0,  the 
boundary  dQ  is  divided  into  three  non-overlapping  parts: 

dn  =rturburv, 


rt  =  {x  G  df2  |n(x)  •  b  <  0}  , 
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^  =  {x  G  dQ  |n(x)  •  b  >  0} ,  Tv  =  {x  G  dQ  |n(x)  •  b  =  0}  . 

With  the  above  partition,  the  boundary  condition  for  Y5  is 

^5  = +  onrbx(0,T), 

Y5  =  -^^+5y5c(x),  onrtx(0,T), 

— q  n  =  0,  on  Tv  x  (0 ,  T7), 

wherein  <5y5  h  (x)  and  Sy5  c  (x)  are  small  scalar  perturbation  functions  that  mimic  the 
uneven  temperature  distribution  on  the  solid  surface.  The  initial  conditions  represent 
a  static  free  surface,  with  liquid  in  the  bottom  region  and  vapor  in  the  top  region.  It  is 
worth  emphasizing  that,  in  contrast  to  existing  boiling  models,  there  is  no  artificial 
manipulation  used  to  serve  as  boiling  onset  in  this  model;  the  initial  liquid  and  vapor 
densities  are  uniform  with  no  perturbations. 


4.1  Two-dimensional  nucleate  boiling 

In  this  example,  we  simulate  boiling  flows  in  a  two-dimensional  rectangular  domain 
12  =  (0, 1)  x  (0,0.5).  The  material  parameters  are  chosen  as  We  =  8.401  x  106, 
y  =  1.333,  Ch°il  =  1.150  x  10~4,  and  CbKoil  =  1.725  x  10“5.  The  initial  conditions 
for  this  problem  are 

po(x)  =  0.3660 

u0(x)  =  0, 
flb(x)  =  0.775. 

The  perturbation  for  the  temperature  on  the  boundary  Sy5  h  (x)  and  Sy5  c  (x)  are  uni¬ 
form  random  distributions  and  satisfy 

Sy5h (x)  e  [—5.0  x  10-2,5.0  x  10-2],  Sy5c(x)  €  [-5.0  x  1(T3,5.0  x  10"3]. 

The  spatial  mesh  consists  of  2048  x  1024  quadratic  NURBS  elements.  The  problem 
is  integrated  up  to  the  final  time  T  =  100.0  with  time  step  fixed  as  At  =  5.0  x  10-4. 
In  Figure  1,  snapshots  of  the  density  are  depicted  at  different  time  steps.  It  can  be 
observed  that  tiny  vapor  bubbles  are  generated  at  discrete  sites  of  the  heated  wall 
surface  during  the  initial  times.  The  increase  of  bubble  size  leads  to  the  increase  of 
buoyancy.  At  about  t  =  18.75,  the  first  three  bubbles  get  detached  from  the  bottom. 
More  bubbles  are  generated  on  the  bottom  surface  subsequently.  Interestingly,  small 
droplets  can  be  observed  at  t  =  62.5  and  t  =  100.0  as  a  result  of  the  breakage  of  the 
liquid  film  when  the  vapor  bubbles  reach  the  free  surface.  There  are  30  bubbles 
formed  in  the  time  interval  of  0  <  t  <  100. 
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Fig.  1  Two-dimensional  nucleate  boiling  simulation:  Density  profiles  at  (a)  t  =  0.0,  (b)  t  =  1.25, 
(c)  t  =  18.75,  (d)  t  =  31.25,  (e)  t  =  62.5,  and  (f)  t  =  100.0. 


4.2  Two-dimensional  film  boiling 

In  the  second  example,  the  same  two-dimensional  problem  considered  in  the  preced¬ 
ing  section  is  simulated  again  with  a  different  parameter  Cb^nl .  Here,  the  parameter 
C^°l1  is  chosen  to  be  4.600  x  10-4,  which  is  four  times  larger  than  that  of  the  pre¬ 
vious  example.  Since  the  fluid  motion  in  this  example  is  slower,  the  simulation  is 
integrated  in  time  up  to  T  =  500.0.  All  the  other  conditions  are  identical  to  those 
of  the  previous  case.  In  Figure  2,  snapshots  of  the  density  at  different  time  steps 
are  depicted.  A  thin  vapor  film  is  gradually  generated  at  the  bottom  during  the  early 
stage  of  the  simulation.  As  time  evolves,  the  interface  becomes  unstable  and  there 
are  vapor  bubbles  formed.  From  t  =  200.0  to  t  =  225.0,  the  first  two  vapor  bubbles 
pinch  off  from  the  vapor  film  and  rise  upward  in  ellipsoidal  shapes.  This  process 
repeats  itself  periodically.  By  final  time  t  =  500.0,  there  are  seven  bubbles  detached 
from  the  vapor  film. 


Isogeometric  Phase-field  Simulation  of  Boiling 


9 


(b) 


Fig.  2  Two-dimensional  film  boiling  simulation:  Density  profiles  at  (a)  t  =  0.0,  (b)  t  =  100.0,  (c) 
t  =  175.0,  (d)  t  =  200.0,  (e)  t  =  225.0,  and  (f)  t  =  500.0.  The  generation  of  the  thin  vapor  film  is 
visible  at  t  =  100.0. 


4.3  Three-dimensional  boiling 

As  the  last  example,  we  simulate  the  Navier-Stokes-Korteweg  equations  in  a  three- 
dimensional  domain  £2  =  (0,1)  x  (0,0.5)  x  (0,0.25).  The  material  properties  are 
chosen  as  We  =  6.533  x  105,  7  =  1.333,  Cftil  =  1.289  x  10~4,  and  CbKoil  =  7.732  x 
10-5.  The  initial  conditions  for  this  three-dimensional  problem  are 

Po(x)  =  0.33565- 0.26675  tanh  ('Q  ~2°''5  , 

u0(x)  =  0, 

y5,0(x)  =  -1.2334 -0.0569  tanh  ('Q  ~2°'15  v7^)  . 

The  perturbations  of  the  temperature  on  the  boundary  8y5  h  (x)  and  Sy5  c  (x)  are  uni¬ 
form  random  distributions  and  satisfy 

Sy5/I(x)  e  [—5.0  x  10_2,5.0 x  10-2],  SySc(x)  e  [—5. Ox  10_3,5.0x  10-3]. 
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The  spatial  mesh  consists  of  600  x  300  x  150  quadratic  NURBS  elements.  The  prob¬ 
lem  is  integrated  in  time  up  to  T  =  20.0  with  a  fixed  time  step  size  At  =  2.0  x  10-3. 
In  Figure  3,  snapshots  of  density  isosurfaces  and  velocity  streamlines  are  presented. 
At  the  initial  stage,  there  is  an  unstable  vapor  film  formed  over  the  heated  wall  sur¬ 
face.  This  film  soon  separates  into  isolated  vapor  bubbles  located  at  random  sites. 
Since  the  simulation  domain  is  very  shallow  in  the  vertical  direction,  these  bubbles 
reach  the  free  surface  before  they  get  fully  detached  from  the  bottom.  When  these 
high-temperature  vapor  bubbles  reach  the  cooled  top  surface,  they  condense  into 
liquid  droplets  instantaneously  (see  Fig.  3  (e)).  At  t  =  20.0,  a  second  round  of  vapor 
bubbles  is  clearly  generated  on  the  bottom  and  the  liquid  droplets  on  the  top  surface 
merge  together. 


Fig.  3  Three-dimensional  boiling  simulation:  Density  isosurfaces  and  velocity  streamlines  at  (a) 
t  =  0.0,  (b)  t  =  0.6,  (c)  t  =  5.0,  (d)  t  =  11.0,  (e)  t  =  14.0,  and  (f)  t  =  20.0. 
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5  Conclusion 

In  this  work,  we  presented  theoretical  and  numerical  methodologies  for  the  study 
of  boiling,  capable  of  describing  complicated  phase  transition  phenomena  without 
resorting  to  empirical  assumptions.  Our  algorithm  is  provably  entropy- stable  and 
second-order  accurate  in  time.  It  provides  a  unified  predictive  framework  for  nu¬ 
cleate  and  film  boiling  in  two  and  three  dimensions.  In  the  future,  the  presented 
methodologies  will  be  applied  to  the  study  of  other  important  phase  transition  phe¬ 
nomena,  such  as  cavitation,  spray  and  mist  formation. 
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